options("digits" =4)

coul <- c("#E0E0E0", "#FFFF99", "#CCFF99", "#00CC00", "#66FFFF","#0066CC", "#006633", "#FF99FF", "#994C00","#FF9999","#FF0000")
coul11 <- c(coul, "#FAE1C7", "#000000", "#FF8000", "#FAA4A4", "#d18975", "#2d543d", "#ff0055", "#5DA5DA", "#B276B2", "#F17CB0", "#009E73", "#F0E442", "#D55E00")  #11 couleurs en ajouter en fct du nm de SM

my_palette = c("D0"= "#E0E0E0",
               "D4"= "#ccFF99",
               "D7"= "#0066CC",
               "D10"= "#006633",
               "D14"= "#FF99FF",
               "D21"= "#994C00",
               "D32"= "#FF9999",
               "D90"= "#FF0000"
               )

theme_CA <- function (base_size = 11, base_family = "", base_line_size = base_size/22, 
                      base_rect_size = base_size/22) 
{
  
  theme_bw(base_size = base_size, base_family = base_family, 
           base_line_size = base_line_size, base_rect_size = base_rect_size) %+replace% 
    theme(panel.border = element_rect(fill = NA, colour = "black", 
                                      linewidth = rel(1)), panel.grid.major = element_blank(),
          plot.title = element_text(size=10), 
          panel.grid.minor = element_blank(), axis.line = element_line(colour = "black", linewidth = rel(1)), legend.key = element_blank(), 
          strip.background = element_rect(fill = "white", colour = "black", 
                                          linewidth = rel(2)), complete = TRUE)
}

seu <- readRDS("/Users/victor/Documents/Master_2/Projet_long/Flux_metabo/results/seu__flux.RDS")

metabolite <- read.csv("/Users/victor/Documents/Master_2/Projet_long/Flux_metabo/data/Results_scFEA/mouse_balance_gut.csv", header = T, row.names = 1)

metabolite <- data.matrix(metabolite)
metabolite0 <- t(metabolite)
seu[['METABOLITES']] <- CreateAssayObject(counts = metabolite0)
DefaultAssay(seu) <- 'METABOLITES'
seu

## An object of class Seurat 
## 29314 features across 30626 samples within 4 assays 
## Active assay: METABOLITES (70 features, 0 variable features)
##  3 other assays present: RNA, SCT, FLUX
##  5 dimensional reductions calculated: pca, umap, pca.flux, umap.metabo, umap.flux

Profils de concentrations

day <- seu@meta.data$orig.ident
metabolite <- as_tibble(metabolite)
metabolite$day <- day
metabolite <- relocate(metabolite, day, .before = AMP)

dt_c0 <- as.matrix(seu@assays$METABOLITES@counts)
yyy <- seu@meta.data$orig.ident 
levels(as.factor(seu@meta.data$orig.ident))

## [1] "D0"  "D4"  "D7"  "D10" "D14" "D21" "D32"

coul_rev <- rev(coul)
for (jj in 1:nrow(dt_c0)) {
    xxx <- dt_c0[jj, ] 
    final_df <- data.frame(var = paste('X', 1:length(xxx), sep=''),  # le même que plus haut
                       meta = xxx,
                       day = yyy)
    order <- rev(c("D0", "D3", "D4", "D5", "D6", "D7", "D10", "D14", "D21", "D32", "D90"))
    final_df <- transform(final_df, day = factor(day, levels = order))
    title <- rownames(dt_c0)[jj] # récupère le noms du module
    aa <- ggplot(final_df, aes(x = meta, y = day, fill = day)) +
      geom_density_ridges() +
      theme_ridges() +
      theme(legend.position = 'none') +
      ggtitle(title) +
      theme(plot.title = element_text(hjust = 0.5)) +
      scale_fill_manual(values = coul_rev) + 
      xlab("Concentrations") +
      geom_vline(xintercept = median(final_df$meta))
    plot(aa)
   }

## Picking joint bandwidth of 1.21e-05

## Picking joint bandwidth of 0.00141

## Picking joint bandwidth of 0.00019

## Picking joint bandwidth of 0.000509

## Picking joint bandwidth of 6.22e-05

## Picking joint bandwidth of 1.46e-07

## Picking joint bandwidth of 5.92e-07

## Picking joint bandwidth of 5.74e-06

## Picking joint bandwidth of 0.00111

## Picking joint bandwidth of 1.47e-06

## Picking joint bandwidth of 2.82e-07

## Picking joint bandwidth of 0.000479

## Picking joint bandwidth of 9.4e-08

## Picking joint bandwidth of 3.7e-06

## Picking joint bandwidth of 0.000644

## Picking joint bandwidth of 0.000598

## Picking joint bandwidth of 3.11e-08

## Picking joint bandwidth of 9.42e-06

## Picking joint bandwidth of 3.55e-05

## Picking joint bandwidth of 0.000146

## Picking joint bandwidth of 0.000848

## Picking joint bandwidth of 5.91e-06

## Picking joint bandwidth of 1.07e-05

## Picking joint bandwidth of 1.53e-05

## Picking joint bandwidth of 0.00156

## Picking joint bandwidth of 3.28e-06

## Picking joint bandwidth of 0.0066

## Picking joint bandwidth of 5.83e-07

## Picking joint bandwidth of 0.00427

## Picking joint bandwidth of 4.94e-06

## Picking joint bandwidth of 4.86e-05

## Picking joint bandwidth of 2.18e-09

## Picking joint bandwidth of 0.000131

## Picking joint bandwidth of 9.73e-06

## Picking joint bandwidth of 5.42e-06

## Picking joint bandwidth of 0.000631

## Picking joint bandwidth of 1.88e-06

## Picking joint bandwidth of 3.74e-05

## Picking joint bandwidth of 4.53e-06

## Picking joint bandwidth of 1.59e-06

## Picking joint bandwidth of 2.16e-09

## Picking joint bandwidth of 4.83e-06

## Picking joint bandwidth of 1.02e-08

## Picking joint bandwidth of 0.00121

## Picking joint bandwidth of 3.62e-05

## Picking joint bandwidth of 8.67e-06

## Picking joint bandwidth of 0.000136

## Picking joint bandwidth of 0.000112

## Picking joint bandwidth of 6.58e-08

## Picking joint bandwidth of 9.84e-06

## Picking joint bandwidth of 2.55e-06

## Picking joint bandwidth of 0.00013

## Picking joint bandwidth of 1.56e-06

## Picking joint bandwidth of 6.86e-07

## Picking joint bandwidth of 0.000152

## Picking joint bandwidth of 0.00165

## Picking joint bandwidth of 7.17e-05

## Picking joint bandwidth of 1.88e-06

## Picking joint bandwidth of 0.00506

## Picking joint bandwidth of 0.000129

## Picking joint bandwidth of 0.000598

## Picking joint bandwidth of 1.56e-08

## Picking joint bandwidth of 2.89e-06

## Picking joint bandwidth of 3.23e-09

## Picking joint bandwidth of 6.1e-07

## Picking joint bandwidth of 2.11e-07

## Picking joint bandwidth of 9e-05

## Picking joint bandwidth of 9e-05

## Picking joint bandwidth of 1.06e-06

## Picking joint bandwidth of 0.000259

Cinétique

means <- metabolite %>%
    group_by(day) %>%
    summarise(across(everything(), mean, na.rm= TRUE)) %>%
    t()

## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(everything(), mean, na.rm = TRUE)`.
## ℹ In group 1: `day = D0`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))

means <- cbind(rownames(means), means)
colnames(means) <- as.character(means[1,])
colnames(means)[1] <- "metabolites"
means <- means[-1,]
means <- as_tibble(means)
means

## # A tibble: 70 × 8
##    metabolites             D0         D4         D7      D10   D14   D21   D32  
##    <chr>                   <chr>      <chr>      <chr>   <chr> <chr> <chr> <chr>
##  1 AMP                     -0.001827  -0.001888  "-0.00… "-0.… -0.0… -0.0… -0.0…
##  2 Pyruvate                0.02118    0.10085    "0.054… "0.0… 0.02… 0.03… 0.02…
##  3 Acetyl.CoA              -0.01212   -0.02067   "-0.01… "-0.… -0.0… -0.0… -0.0…
##  4 Glutamate               -0.002500  -0.003629  " 0.00… "-0.… -0.0… -0.0… -0.0…
##  5 X2OG                    0.007885   0.010517   "0.008… "0.0… 0.00… 0.00… 0.00…
##  6 Oxaloacetate            -0.0005159 -0.0005223 "-0.00… "-0.… -0.0… -0.0… -0.0…
##  7 Glycine                 0.005076   0.005048   "0.005… "0.0… 0.00… 0.00… 0.00…
##  8 Succinate               0.002201   0.002276   "0.002… "0.0… 0.00… 0.00… 0.00…
##  9 UDP.N.acetylglucosamine -0.005280  -0.001113  " 0.02… " 0.… -0.0… -0.0… -0.0…
## 10 lysine                  0.01197    0.01194    "0.011… "0.0… 0.01… 0.01… 0.01…
## # ℹ 60 more rows

maxi <-  c()
for (i in 1:70) {
  maxi <- c(maxi, max(as.numeric(means[i, 2:8])))
}
mini <-  c()
for (i in 1:70) {
  mini <- c(mini, min(as.numeric(means[i, 2:8])))
}
delta <- maxi - mini
means$conc_min <- mini
means$conc_max <- maxi
means$conc_delta <- delta
means

## # A tibble: 70 × 11
##    metabolites       D0    D4    D7    D10   D14   D21   D32   conc_min conc_max
##    <chr>             <chr> <chr> <chr> <chr> <chr> <chr> <chr>    <dbl>    <dbl>
##  1 AMP               -0.0… -0.0… "-0.… "-0.… -0.0… -0.0… -0.0… -1.93e-3 -1.82e-3
##  2 Pyruvate          0.02… 0.10… "0.0… "0.0… 0.02… 0.03… 0.02…  2.12e-2  1.01e-1
##  3 Acetyl.CoA        -0.0… -0.0… "-0.… "-0.… -0.0… -0.0… -0.0… -2.07e-2 -1.21e-2
##  4 Glutamate         -0.0… -0.0… " 0.… "-0.… -0.0… -0.0… -0.0… -4.55e-3  2.28e-3
##  5 X2OG              0.00… 0.01… "0.0… "0.0… 0.00… 0.00… 0.00…  7.89e-3  1.05e-2
##  6 Oxaloacetate      -0.0… -0.0… "-0.… "-0.… -0.0… -0.0… -0.0… -5.22e-4 -5.16e-4
##  7 Glycine           0.00… 0.00… "0.0… "0.0… 0.00… 0.00… 0.00…  5.05e-3  5.08e-3
##  8 Succinate         0.00… 0.00… "0.0… "0.0… 0.00… 0.00… 0.00…  2.20e-3  2.28e-3
##  9 UDP.N.acetylgluc… -0.0… -0.0… " 0.… " 0.… -0.0… -0.0… -0.0… -5.28e-3  2.74e-2
## 10 lysine            0.01… 0.01… "0.0… "0.0… 0.01… 0.01… 0.01…  1.19e-2  1.20e-2
## # ℹ 60 more rows
## # ℹ 1 more variable: conc_delta <dbl>

means <- arrange(means, desc(conc_delta))
ggplot(means, aes(conc_delta)) + geom_dotplot() + scale_x_log10()

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

means_2 <- filter(means, conc_delta > 1e-3) ## Reste 32 métabolites

Profils de concentrations finaux

T_F <- rownames(dt_c0) %in% means_2$metabolites
dt <- dt_c0[T_F == T,]
coul_rev <- rev(coul)
for (jj in 1:nrow(dt)) {
    xxx <- dt[jj, ] 
    final_df <- data.frame(var = paste('X', 1:length(xxx), sep=''),  # le même que plus haut
                       meta = xxx,
                       day = yyy)
    order <- rev(c("D0", "D3", "D4", "D5", "D6", "D7", "D10", "D14", "D21", "D32"))
    final_df <- transform(final_df, day = factor(day, levels = order))
    title <- rownames(dt)[jj] # récupère le noms du module
    aa <- ggplot(final_df, aes(x = meta, y = day, fill = day)) +
      geom_density_ridges() +
      theme_ridges() +
      theme(legend.position = 'none') +
      ggtitle(title) +
      theme(plot.title = element_text(hjust = 0.5)) +
      scale_fill_manual(values = coul_rev) + 
      xlab("Concentrations") +
      geom_vline(xintercept = median(final_df$meta))
    plot(aa)
   }

## Picking joint bandwidth of 0.00141

## Picking joint bandwidth of 0.00019

## Picking joint bandwidth of 0.000509

## Picking joint bandwidth of 6.22e-05

## Picking joint bandwidth of 0.00111

## Picking joint bandwidth of 0.000479

## Picking joint bandwidth of 0.000644

## Picking joint bandwidth of 0.000598

## Picking joint bandwidth of 3.55e-05

## Picking joint bandwidth of 0.000146

## Picking joint bandwidth of 0.000848

## Picking joint bandwidth of 0.00156

## Picking joint bandwidth of 0.0066

## Picking joint bandwidth of 0.00427

## Picking joint bandwidth of 4.86e-05

## Picking joint bandwidth of 0.000131

## Picking joint bandwidth of 0.000631

## Picking joint bandwidth of 3.74e-05

## Picking joint bandwidth of 0.00121

## Picking joint bandwidth of 3.62e-05

## Picking joint bandwidth of 0.000136

## Picking joint bandwidth of 0.000112

## Picking joint bandwidth of 0.00013

## Picking joint bandwidth of 0.000152

## Picking joint bandwidth of 0.00165

## Picking joint bandwidth of 7.17e-05

## Picking joint bandwidth of 0.00506

## Picking joint bandwidth of 0.000129

## Picking joint bandwidth of 0.000598

## Picking joint bandwidth of 9e-05

## Picking joint bandwidth of 9e-05

## Picking joint bandwidth of 0.000259

Clustering

seu <- FindVariableFeatures(seu, selection.method = "vst", verbose = F)

## Warning in eval(predvars, data, env): NaNs produced

## Warning in hvf.info$variance.expected[not.const] <- 10^fit$fitted: number of
## items to replace is not a multiple of replacement length

seu

## An object of class Seurat 
## 29314 features across 30626 samples within 4 assays 
## Active assay: METABOLITES (70 features, 70 variable features)
##  3 other assays present: RNA, SCT, FLUX
##  5 dimensional reductions calculated: pca, umap, pca.flux, umap.metabo, umap.flux

seu <- ScaleData(seu, features = rownames(seu), assay = 'METABOLITES', verbose = F)
seu <- RunPCA(seu, features = VariableFeatures(object = seu), reduction.name = 'pca.metabo', verbose = F)

## Warning in irlba(A = t(x = object), nv = npcs, ...): You're computing too large
## a percentage of total singular values, use a standard svd instead.

## Warning: Cannot add objects with duplicate keys (offending key: PC_), setting
## key to 'pca.metabo_'

DimPlot(seu, dims = 1:2,reduction = "pca.metabo", group.by = "orig.ident", label = T)

ElbowPlot(seu, ndims = 30, reduction = 'pca.metabo')   # 9 dim

Clustree and clusterization

seu <- RunUMAP(seu, dims = 1:9, assay = 'METABOLITES', reduction = 'pca.metabo', reduction.name = "umap.metabo", verbose = F)

## Warning: Cannot add objects with duplicate keys (offending key: UMAP_) setting
## key to original value 'umap.metabo_'

seu <- FindNeighbors(seu, dims = 1:9, reduction = 'pca.metabo', verbose = F, assay = 'METABOLITES')
seu <- FindClusters(seu , graph.name = 'METABOLITES_snn', verbose = F)

DimPlot(seu, reduction = "umap.metabo", label = T, label.box = T, group.by = "METABOLITES_snn_res.0.3")

DimPlot(seu, reduction = "umap.metabo", label = T, label.box = T, group.by = "orig.ident")

DimPlot(seu, reduction = "umap", label = T, label.box = T, group.by = "METABOLITES_snn_res.0.3")

Barplot

cell <- cbind(Cells(seu), seu@meta.data$orig.ident)
colnames(cell) <- c("barcode", "dataset")
day <- seu@meta.data$orig.ident
cell <- as_tibble(cell)
cell$day <- day
cell$rna_cluster <- seu@meta.data$SCT_snn_res.0.8
cell$flux_cluster <- seu@meta.data$FLUX_snn_res.0.5
cell$metabo_cluster <- seu@meta.data$METABOLITES_snn_res.0.3

unique(cell$flux_cluster)

##  [1] 12 16 4  11 5  14 10 17 6  3  8  13 7  2  1  15 9  0 
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

unique(cell$rna_cluster)

##  [1] 6  0  14 11 4  9  15 1  8  19 17 18 7  16 3  12 2  5  10 13
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

unique(cell$metabo_cluster)

##  [1] 3  9  6  12 4  7  2  11 0  10 1  8  5 
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12

ggplot(cell, aes(x = day, fill = rna_cluster)) +
  geom_bar() +
  labs(title = "Repartition des cluster RNA en fonction des jours", 
       y = "RNA clusters") +
  theme_CA() +
  scale_fill_manual(values = coul11) +
  theme(axis.title.x = element_blank())

ggplot(cell, aes(x = rna_cluster, fill = day)) +
  geom_bar() +
  labs(title = "Repartition des jours dans les clusters RNA", 
       y = "Jour") +
  theme_CA() +
  scale_fill_manual(values = my_palette) +
  theme(axis.title.x = element_blank())

ggplot(cell, aes(x = day, fill = flux_cluster)) +
  geom_bar() +
  labs(title = "Repartition des cluster de flux en fonction des jours", 
       y = "Flux clusters") +
  theme_CA() +  
  scale_fill_manual(values = coul11) +
  theme(axis.title.x = element_blank())

ggplot(cell, aes(x = flux_cluster, fill = day)) +
  geom_bar() +
  labs(title = "Repartition des jours dans les clusters de flux", 
       y = "Jour") +
  theme_CA() +
  scale_fill_manual(values = my_palette) +
  theme(axis.title.x = element_blank())

ggplot(cell, aes(x = day, fill = metabo_cluster)) +
  geom_bar() +
  labs(title = "Repartition des cluster de flux en fonction des jours", 
       y = "Flux clusters") +
  theme_CA() +  
  # scale_fill_manual(values = coul11) +
  theme(axis.title.x = element_blank())

ggplot(cell, aes(x = metabo_cluster, fill = day)) +
  geom_bar() +
  labs(title = "Repartition des jours dans les clusters de métabolite", 
       y = "Jour") +
  theme_CA() +
  scale_fill_manual(values = my_palette) +
  theme(axis.title.x = element_blank())

Analyse_metabolites

Malassigne_Victor

4/19/2023

Profils de concentrations

Cinétique

Profils de concentrations finaux

Clustering

Clustree and clusterization

Barplot